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Abstract.  The  ground  state  of  the  system  of  charged  particles  of  one  species  confined  by  the  three- 
dimensional,  isotropic,  and  parabolic  potential  is  investigated  by  molecular  dynamics  simulations. 
It  is  shown  that,  with  the  increase  of  the  system  size  or  the  number  of  particles  in  the  system  N, 
the  ground  state  changes  from  the  shell-structured  system  to  the  finite  bcc  lattice  with  reconstructed 
surface.  The  critical  value  of  the  transition  is  estimated  to  be  between  N = 104  and  N = 2 x 104. 
The  nucleation  of  the  bcc  lattice  in  the  shell-structured  cluster  of  2 x 104  ions  is  observed. 


INTRODUCTION 

The  systems  of  charges  of  one  species  confined  in  the  external  field  provide  us  with 
a stage  where  various  properties  of  strongly  coupled  plasmas  manifest  themselves  in 
a simple  and  clear  way.  This  system  can  be  realized  in  the  Penning-Malmberg  [l]and 
Paul  traps  [2,  3]  and  direct  experimental  observations  of  these  properties  have  been  an 
excitement  in  recent  years. 

The  ground  state  of  these  finite  systems  is  strongly  influenced  by  the  geometry  of 
confinement.  When  the  external  potential  is  the  spherically  symmetric  parabola  and  the 
number  of  charges  in  the  system  is  small,  the  structures  such  as  icosahedron  or  spherical 
shells  give  the  ground  state.  When  the  system  has  sufficiently  large  number  of  particles, 
it  is  expected  that  the  ground  state  contains  the  bcc  lattice  as  its  main  part:  The  ground 
state  of  the  (infinite)  one-component  plasma  is  the  bcc  structure. 

In  experiments,  the  bcc  lattice  at  the  central  part  has  been  observed  for  large  cloud  of 
charges  of  the  size  more  than  105  [4,  5J.  The  transition  of  the  ground  state  from  the  shell 
structure  to  the  bcc  lattice  has  been  expected  to  occur  at  the  system  size  of  2 x 105  or 
larger  [6],  It  seems,  however,  that  numerical  simulations  of  larger  systems  may  give  more 
precise  information  related  to  this  transition.  The  purpose  of  this  paper  is  to  compare  the 
cohesive  energy,  sometimes  called  the  Madelung  energy,  of  these  structures  and  show 
that  the  transition  occurs  between  N = 104  and  N = 2 x 104  by  molecular  dynamics 
simulations  of  large  clusters. 
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PARAMETERS  AND  MADELUNG  ENERGY 

Let  us  here  fix  some  notations  for  these  clusters  of  charges.  The  Hamiltonian  of  our 
system  II  is  given  by  H = K + U,  where  K is  the  kinetic  energy 


A I / dr, 
A = > -?/?  — - 
St  2 \dt 


(I) 


and  U,  the  potential  energy 
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We  define  the  Wigner-Seitz  radius  «it\s-  by 
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and  rewrite  K and  U into  the  form 
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Here  the  coordinates  of  particles  are  normalized  by  «ir.s  and  the  time  is  normalized  by 
to. 

In  the  fluid  approximation  where  charges  are  regarded  as  continuum,  the  ground  state 
is  the  uniform  distribution  up  to  the  radius  /? 
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with  the  density 
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We  note  that  uh-.s  defined  by  (3)  is  the  average  mean  distance  or  the  ion-sphere  radius 
in  this  approximation; 
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The  potential  energy  of  the  ground  state  in  the  fluid  limit  is  given  by 
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For  given  distribution  of  particles,  we  define  the  Madelung  energy  by 

U M = f ^ homo  ■ 


(10) 


At  zero  temperature,  the  behavior  of  the  system  is  determined  only  by  U.  Therefore 
(5)  indicates  that  our  system  is  characterized  by  only  one  parameter,  the  system  size  N. 
At  finite  temperatures,  we  define  the  parameter  T which  represents  the  strength  of  the 
Coulomb  coupling  by 


r = 
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The  static  properties  of  the  system  at  finite  temperatures  are  characterized  by  the  param- 
eter T and  the  system  size  N. 

In  the  fluid  approximation,  the  radius  of  uniform  distribution  oscillates  with  the 
frequency 
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This  frequency  characterizes  the  macroscopic  evolution  of  our  system.  One  of  character- 
istic time  scales  of  the  microscopic  evolution  may  be  the  time  to  move  the  mean  distance 
by  thermal  velocity  aWsl{kBT/m)l/ 2 which  is  related  to  ujp  as 
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NUMERICAL  METHOD 
Scaling  with  system  size 

We  perform  the  molecular  dynamics  simulations.  In  this  case,  most  of  the  computa- 
tional time  is  consumed  in  the  computation  of  the  force  on  each  particle  and,  with  the 
naive  method,  the  computational  time  increases  rapidly  in  proportion  to  N2  in  our  sys- 
tem where  particles  interact  via  the  long-ranged  Coulomb  potential.  We  therefore  adopt 
the  fast  multipole  method  [7]  which  enables  one  to  perform  molecular  dynamics  simu- 
lations of  long-range  force  systems  with  the  O(N)- scaling  under  controlled  accuracy. 

In  the  fast  multipole  method,  the  system  is  divided  recursively  into  small  cells  and  the 
interaction  between  particles  belonging  to  well-separated  cells  is  computed  based  on  the 
multipole  expansion  and  Taylor  expansion.  The  number  of  smallest  cells  is  8level  and  the 
level  is  adjusted  so  that  the  smallest  cell  has  about  8 particles  on  the  average  [7].  For 
the  system  size  of  N = 10s,  we  adopt  the  following  parameters  for  the  fast  multipole 
method:  the  level  of  5,  multipole  expansion  up  to  26-th  order,  the  Taylor  expansion 
up  to  6th  order,  the  well-separatedness  of  2 (in  the  highest  level,  direct  computation  is 
applied  for  nearest  and  next-nearest  neighbor  cells).  In  the  course  of  molecular  dynamics 
simulation,  the  accuracy  of  the  Madelung  energy  is  kept  at  least  for  the  first  4 digits  and 
the  direct  computation  is  performed  for  final  relaxed  states.  The  parameters  for  smaller 
systems  are  also  determined  to  give  the  accuracy  at  least  of  the  same  order. 
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FIGURE  I.  Madelung  energy  vs.  system  size  for  spherical  hoc  matter.  Solid  circles  are  the  values  with 
relaxed  surface  and  solid  line  is  the  interpolation.  Open  and  gray  circles  arc  examples  of  values  without 
relaxation. 


FIGURE  2.  Particle  distribution  near  the  equatorial  plane  in  spherical  bcc  matter  with  relaxed  surface 
for  N = 120.032.  Positions  of  particles  in  the  domain  |c|  < 2.7«n\s-  are  projected. 
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The  temperature  is  controlled  by  the  Nose-Hoover  method.  In  order  to  keep  the 
homogeneity  of  the  temperature,  we  attach  multiple  thermostats  each  controlling  the 
kinetic  energy  of  about  5000  particles. 


Initial  conditions 

For  the  system  of  the  size  N,  the  radius  of  particle  distribution  in  the  fluid-limit  ground 
state  is  determined  as  (6).  The  initial  positions  for  the  shell  structure  is  given  by  the 
uniform  random  distribution  of  particles  within  this  sphere.  The  velocities  are  given 
by  the  random  distribution  corresponding  to  the  temperature  specified  by  the  value  of 
T.  After  annealing  with  T values  between  100  and  150  for  sufficiently  many  steps  in 
measures  of  the  time  scales  (12)  and  (13),  we  slowly  decrease  the  temperature. 

Another  set  of  initial  conditions  is  the  spherical  cutout  of  the  three-dimensional  lattice 
and  the  random  distribution  of  velocities.  In  this  case,  the  value  of  T is  kept  just  above 
the  melting  point  (low  temperature  side)  of  the  infinite  lattice  for  a sufficiently  long  time 
to  allow  the  system  to  relax  before  slowly  increased. 


MADELUNG  ENERGY  OF  RELAXED  SPHERICAL  BCC 

MATTER 

The  spherical  cutout  of  he  bcc  lattice  has  rather  high  surface  energy  as  is  shown  in  Fig.  1 . 
When  we  anneal  this  system  keeping  the  temperature  near  but  below  the  melting  point, 
we  obtain  the  spherical  bcc  matter  with  reconstructed  surface.  The  Madelung  energy  of 
the  spherical  bcc  matter  with  reconstructed  surface  is  also  shown  in  Fig.  I for  the  system 
sizes  N = 4544,  10,464,  20,288,  48,928,  and  120,032.  We  observe  that  the  relaxation 
near  the  surface  largely  reduces  the  Madelung  energy  [8].  This  is  not  expected  from 
the  high  values  of  the  Madelung  energy  of  unrelaxed  spherical  bcc  matter.  An  example 
of  the  structure  is  shown  in  Figs.  2 and  3.  We  observe  that  the  relaxation  has  occurred 
within  a few  layers  at  the  surface 

In  the  limit  of  very  large  values  of  N,  the  Madelung  energy  may  have  the  form 

Um  Es 

Nq2/aWS  = °°+NW'  (14) 

where  the  second  term  expresses  the  effect  of  the  surface.  When  fitted  to  this  form  with 
Eao  = -0.895929  (the  Madelung  energy  of  an  infinite  bcc  solid),  as  shown  in  Fig.  I by 
the  solid  line,  we  have  Es  = 0.0G05. 


MADELUNG  ENERGY  OF  SHELL-STRUCTURED  SPHERICAL 

CLUSTERS 

Starting  from  the  random  distribution  in  the  sphere  of  radius  R given  by  (6),  we  first  an- 
neal the  system  at  the  temperature  above  the  melting  temperature  of  bulk  one-component 
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FIGURE  5.  Particle  distribution  in  shell-structured  spherical  cluster  with  N = 100,000. 


plasma.  The  typical  value  of  T is  around  100.  We  then  slowly  decrease  the  temperature 
of  the  system.  The  results  are  shown  in  Figs.  4,  5,  and  6. 

The  values  of  the  Madelung  energy  of  the  shell-structured  spherical  system  of  the  size 
Ar  = 5x  10\  101,  20288,  and  10r>  are  shown  in  Fig.  7.  These  values  of  the  Medelung 
energy  can  be  interpolated  by  the  formula 

— ~ — = —0.8910  + 0.0219AT1/3.  (15) 
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The  largest  simulated  system  in  reports  was  N = 2 x 10*  [9], 


TRANSITION  OF  THE  GROUND  STATE 

The  Madelung  energies  of  the  shell-structured  clusters  and  of  finite  bcc  matters  with 
relaxed  surface  are  compared  in  Fig.  7.  We  observe  that,  when  the  system  size  exceeds 
the  critical  value  Nc,  the  ground  state  changes  from  the  shell-structured  cluster  to  the 
finite  bcc  lattice  with  relaxed  surface.  The  critical  system  sized  is  given  by  the  crossing  of 
these  Madelung  energies  and,  in  Fig.  7,  is  estimated  to  be  around  101.  We  here  note  that 
these  Madelung  energy,  especially  of  shell  structures,  could  become  lower  by  annealing 
the  system  repeatedly.  The  crossing,  however,  does  not  seem  to  be  larger  than  2 x 101. 
Taking  into  account  this  possibility,  we  here  conclude  that 

101  < Nc  < 2x  101.  (16) 


NUCLEATION  OF  BCC  LATTICE 

The  shell-structured  cluster  is  still  is  a local  minimum  of  the  Madelung  energy  in  the 
domain  where  N > Nc  and  the  bcc  lattice  with  relaxed  surface  is  the  global  minimum. 
Therefore  it  is  not  strange  that  the  shell-structured  cluster  does  not  evolve  to  the  bcc 
lattice  with  relaxed  surface  in  simulations.  One  may,  however,  expect  that  there  is  a 
chance  for  the  bcc  lattice  develops  in  the  clusters  of  larger  sizes.  We  here  show  one  of 
such  examples. 

The  distribution  of  particles  near  the  equatorial  plane  of  the  cluster  of  N = 20,288 
is  shown  in  the  left  panel  of  Fig.  8.  We  observe  that  there  exist  a part  where  we  have 
the  structure  with  the  straight  planes  rather  than  a curved  plane  in  accordance  with  the 
surface.  The  structure  factor  for  the  particles  in  this  domain  is  shown  in  the  right  panel 
of  Fig.  8.  We  clearly  see  that  the  Bragg  spots  forms  the  fee  structure  in  the  wave-number 
space.  We  have  confirmed  that  the  spacing  between  Bragg  spots  agree  with  those  of  the 
bcc  lattice  with  the  average  density  of  this  cluster.  This  domain  thus  forms  the  bcc  lattice 
nucleated  from  the  shell-structured  cluster. 
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FIGURE  7.  Madelung  energy  of  spherical  coulomb  system  vs.  N 93  Filled  circles  are  relaxed  spher 
ical  bcc  matter  and  squares  are  shell  structured  systems. 


* * 


FIGURE  8.  Bragg  pattern  from  the  central  part  (circled)  of  the  shell-structured  cluster  of  N = 20, 288. 
Shown  are  consecutive  three  planes  in  reciprocal  space  where  spots  form  the  fee  lattice  corresponding  to 
bcc  structure  in  real  space. 
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CONCLUSIONS 


It  is  shown  that  the  ground  state  of  the  spherical  coulomb  clusters  changes  from  shell- 
structured  clusters  to  the  finite  bee  lattice  with  relaxed  surface  at  the  system  size  between 
104  and  2 x If)1,  In  the  shell-structured  cluster  of  2 x If)4,  the  nucleation  of  the  bcc  lattice 
is  observed. 
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